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Scientific  Progress 


We  have  performed  research  on  the  design  of  new  algorithms  and  improvement  of  existing  algorithms  for  high  order  accurate 
finite  difference  and  finite  volume  weighted  essentially  non-oscillatory  (WENO)  schemes  and  discontinuous  Galerkin  finite 
element  methods,  for  solving  partial  differential  equations  with  discontinuous,  rapidly  changing,  or  multiscale  solutions. 
Applications  to  computational  fluid  dynamics,  astrophysical  problems,  and  pedestrian  flows  are  addressed.  The  objective  of 
improving  the  range  of  applicability,  efficiency,  robustness,  and  scalability  in  massive  parallel  environment  of  the  proposed 
methods  for  various  physical  problems  has  been  achieved.  Particular  attention  has  been  paid  to  army  related  applications 
including  the  pedestrian  flow  problems.  There  are  18  refereed  journal  papers  published  which  have  quoted  partial  support  from 
this  ARO  grant. 

We  explicitly  construct  the  entropy  solutions  for  the  Lighthill-Whitham-Richards  (LWR)  traffic  flow  model  with  a  flow-density 
relationship  which  is  piecewise  quadratic,  concave,  but  not  continuous  at  the  junction  points  where  two  quadratic  polynomials 
meet,  and  with  piecewise  linear  initial  condition  and  piecewise  constant  boundary  conditions.  The  existence  and  uniqueness  of 
entropy  solutions  for  such  conservation  laws  with  discontinuous  fluxes  are  not  known  mathematically.  We  have  used  the 
approach  of  explicitly  constructing  the  entropy  solutions  to  a  sequence  of  approximate  problems  in  which  the  flow-density 
relationship  is  continuous  but  tends  to  the  discontinuous  flux  when  a  small  parameter  in  this  sequence  tends  to  zero.  The  limit 
of  the  entropy  solutions  for  this  sequence  is  explicitly  constructed  and  is  considered  to  be  the  entropy  solution  associated  with 
the  discontinuous  flux.  We  apply  this  entropy  solution  construction  procedure  to  solve  four  representative  traffic  flow  cases, 
compare  them  with  numerical  solutions  obtained  by  a  high  order  weighted  essentially  non-oscillatory  (WENO)  scheme,  and 
discuss  the  results  from  traffic  flow  perspectives. 

High  order  accurate  weighted  essentially  non-oscillatory  (WENO)  schemes  have  been  used  extensively  in  numerical  solutions 
of  hyperbolic  partial  differential  equations  and  other  convection  dominated  problems.  However  the  WENO  procedure  can  not  be 
applied  directly  to  obtain  a  stable  scheme  when  negative  linear 
weights  are  present.  We  have  performed  a  detailed  study  on  the  positivity  of 

linear  weights  in  a  few  typical  WENO  procedures,  including  WENO  interpolation,  WENO  reconstruction  and  WENO 
approximation  to  first  and  second  derivatives,  and  WENO  integration.  Explicit  formulae  for  the  linear  weights  are  also  given  for 
these  WENO  procedures.  These  results  should  be  useful  for  future  design  of  WENO  schemes  involving  interpolation, 
reconstruction,  approximation  to  first  and  second  derivatives,  and  integration  procedures. 

The  interaction  between  a  shock  wave  and  two  counter  rotating  vortices  is  simulated  systematically  through  solving  the  two 
dimensional,  unsteady  compressible  Navier-Stokes  equations  using  a  fifth  order  weighted  essentially  nonoscillatory  (WENO) 
finite  difference  scheme.  The  main  purpose  of  this  study  is  to  reveal  the  mechanism  of  sound  generation  in  the  interaction 
between  a  shock  wave  and  two  counter  rotating  vortices.  It  is  found  that  there  are  two  regimes  of  sound  generation  in  this 
interaction.  The  first  regime  corresponds  to  the  shock  interaction  with  two  isolated  vortices,  in  which  the  sound  wave  generated 
by  the  interaction  between  the  shock  wave  and  two  counter  rotating 

vortices  equals  to  the  linear  combination  of  the  sound  waves  generated  by  the  interactions  between  the  same  shock  wave  and 
each  vortex.  The  second  regime  corresponds  to  the  shock  interaction  with  a  coupled  vortex  pair,  in  which  the  sound  wave 
comes  from  two  processes.  One  is  the  vortex  coupling,  and  the  second  is  the  interaction  between  the  shock  wave  and  the 
coupled  vortex  pair. 

We  have  developed  a  macroscopic  model  for  pedestrian  flow  using  the  dynamic  continuum  modeling  approach.  We  consider  a 
two-dimensional  walking  facility  that  is  represented  as  a  continuum  within  which  pedestrians  can  move  freely  in  any  direction.  A 
pedestrian  chooses  a  route  based  on  his  or  her  memory  of  the  shortest  path  to  the  desired  destination  when  the  facility  is 
empty,  and  at  the 

same  time  tries  to  avoid  high  densities.  In  this  model,  pedestrian  flow  is  governed  by  a  two-dimensional  conservation  law,  and  a 
general  speed-flow-density  relationship  is  considered.  The  model  equation  is  solved  numerically  using  the  discontinuous 
Galerkin  method,  and  a  numerical  example  is  employed  to  demonstrate  both  the  model  and  the  effectiveness  of  the  numerical 
method. 

We  have  presented  further  development  of  the  local  discontinuous  Galerkin  (LDG) 

method  designed  earlier  by  us  and  a  new  dissipative  discontinuous  Galerkin  (DG)  method  for  the  Hunter-Saxton  equation.  The 
numerical  fluxes  for  the  LDG  and  DG  methods  are  based  on  the  upwinding  principle.  The  resulting  schemes  provide  additional 
energy  dissipation  and  better  control  of  numerical  oscillations  near  derivative  singularities.  Stability  and  convergence  of  the 
schemes  are  proved  theoretically,  and  numerical  simulation  results  are  provided  to  compare  with  the  earlier  schemes. 

Using  a  combination  of  critical  point  theory  of  ordinary  differential  equations  and  numerical  simulation  for  the  three-dimensional 
unsteady  Navier-Stokes  equations,  we  study  possible  flow  structures  of  the  vortical  flow,  especially  the  unsteady  vortex 
breakdown  in  the  interaction  between  a  normal  shockwave  and  a  longitudinal  vortex.  The  topological  structure  contains  two 
parts.  One  is  the  sectional  streamline  pattern  in  the  cross-section  perpendicular  to  the  vortex  axis.  The  other  is  the  sectional 
streamline  pattern  in  the  symmetrical  plane.  In  the  cross-section  perpendicular  to  the  vortex  axis,  the  sectional  streamlines 
have  spiral  or  center  patterns  depending  on  a  specific  function.  Depending  on  the  sign  of  a  parameter  in  this  function,  the 


sectional  streamlines  spiral  either  inward  or  outward  in  the  near  region  of  the  center,  or  form  a  nonlinear  center.  If  this 
parameter  changes  its  sign  along  the  vortex  axis,  one  or  more  limit  cycles  appear  in  the  sectional  streamlines  in  the 
cross-section  perpendicular  to  the  vortex  axis.  Numerical  simulation  for  two  typical  cases  of  shock  induced  vortex  breakdown  is 
performed.  The  onset  and  time  evolution  of  the  vortex  breakdown  are  studied.  It  is  found  that  there  are  more  limit  cycles  for  the 
sectional  streamlines  in  the  cross-section  perpendicular  to  the  vortex  axis.  In  addition,  we  find  that  there  are  quadru-helix 
structures  in  the  tail  of  the  vortex  breakdown. 

The  21  cm  emission  and  absorption  from  gaseous  halos  around  the  first  generation  of  stars  substantially  depend  on  the 
Wouthuysen-Field  (W-F)  coupling,  which  relates  the  spin  temperature  with  the  kinetic  temperature  of  hydrogen  gas  via  the 
resonant  scattering  between  Ly$\alpha$  photons  and  neutral  hydrogen.  Therefore,  the  existence  of  Ly$\alpha$  photons  in  the 
21  cm  region  is  essential.  Although  the  center  object  generally  is  a  strong  source  of  Ly$\alpha$  photons,  the  transfer  of 
Ly$\alpha$  photons  in  the  21  cm  region  is  very  inefficient,  as  the  optical  depth  of  Ly$\alpha$  photons  is  very  large. 
Consequently,  the  Ly$\alpha$  photons  $\nu_0$  from  the  source  may  not  be  able  to  transfer  to  the  entire  21  cm  region  timely  to 
provide  the  W-F  coupling.  This  problem  is  especially  important  considering  that  the  lifetime  of  first  stars  generally  is  short.  We 
investigate  this  problem  with  the  numerical  solution  of  the  integrodifferential  equation,  which  describes  the  kinetics  of 
Ly$\alpha$  resonant  photons  in  both  physical  and  frequency  spaces.  We  show  that  the  photon  transfer  process  in  the  physical 
space  is  actually  coupled  to  that  in  the 

frequency  space.  First,  the  diffusion  in  the  frequency  space  provides  a  shortcut  for  the  diffusion  in  the  physical  space.  It  makes 
the  mean  time  for  the  escape  of  resonant  photon  in  optical  depth  $\tau$  media  roughly  proportional  to  the  optical  depth  $\tau$, 
not  $\tauA2$.  Second  and  more  importantly,  the  resonant  scattering  is  effective  in  bouncing  photons  with  frequency  $\nu\neq 
\nu_0$  back  to  $\nu_0$.  This  process  can  quickly  restore  $\nu_0$  photons  and  establish  the  local  Boltzmann  distribution  of  the 
photon  spectrum  around  $\nu_0$.  Therefore,  the  mechanism  of  "escape  via  shortcut"  plus  "bounce  back"  enables  the  W-F 
coupling  to  be  properly  realized  in  the  21  cm  region  around  first  stars.  This  mechanism  also  works  for  photons  injected  into  the 
21  cm  region  by  redshift. 

For  high  order  shock  capturing  numerical  methods,  we  have  introduced  a  new  technique  to  work  with  the  hierarchical 
reconstruction  (HR)  idea  of  Liu,  Shu,  Tadmor  and  Zhang,  to  effectively  reduce  spurious  oscillations  without  local  characteristic 
decomposition  for  numerical  capturing  of  discontinuous  solutions.  The  new  technique  uses  HR  on  partial  neighboring  cells, 
which  lowers  the  order  of  the  remainder  while  maintaining  the  theoretical  order  of  accuracy,  essentially  eliminates 
overshoots/undershoots  for  the  fourth  and  fifth  order  cases  (in  one  dimensional  numerical  examples)  and  reduces  the 
numerical  cost. 

We  have  developed  a  local  discontinuous  Galerkin  (LDG)  method  for  the  generalized  Zakharov  system.  Two  energy 
conservations  of  the  LDG  scheme  are  proved  for  the  generalized  Zakharov  system.  Numerical  experiments  for  the  Zakharov 
system  are  presented  to  illustrate  the  accuracy  and  capability  of  the  methods,  including  accuracy  tests,  plane  waves, 
soliton-soliton  collisions  of  the  standard  and  generalized  Zakharov  system  and  a  two-dimensional  problem. 

We  have  considered  two  commonly  used  classes  of  finite  volume  weighted  essentially  non-oscillatory  (WENO)  schemes  in  two 
dimensional  Cartesian  meshes.  We  compare  them  in  terms  of  accuracy,  performance  for  smooth  and  shocked  solutions,  and 
efficiency  in  CPU  timing.  For  linear  systems  both  schemes  are  high  order  accurate,  however  for  nonlinear  systems,  analysis 
and  numerical  simulation  results  verify  that  one  of  them  (Class  A)  is  only  second 

order  accurate,  while  the  other  (Class  B)  is  high  order  accurate.  The  WENO  scheme  in  Class  A  is  easier  to  implement  and 
costs  less  than  that  in  Class  B.  Numerical  experiments  indicate  that  the  resolution  for  shocked  problems  is 
often  comparable  for  schemes  in  both  classes  for  the  same  building  blocks  and  meshes,  despite  of  the  difference  in  their  formal 
order  of  accuracy.  The  results  in  this  work  may  give  some  guidance  in  the  application  of  high  order  finite  volume  schemes  for 
simulating  shocked  flows. 

With  a  state-of-the-art  numerical  method  for  solving  the  integral-differential  equation  of  radiative  transfer,  we  investigate  the  flux 
of  the  Ly$\alpha$  photon  $\nu_0$  emergent  from  an  optically  thick  halo  containing  a  central  light  source.  Our  focus  is  on  the 
time-dependent  effects  of  the  resonant  scattering.  We  first  show  that  the  frequency  distribution  of  photons  in  the  halo  are 
quickly 

approaching  to  a  locally  thermalized  state  around  the  resonant  frequency,  even  when  the  mean  intensity  of  the  radiation  is 
highly  time-dependent. 

We  have  continued  our  earlier  effort  in  developing  local  discontinuous  Galerkin  (LDG)  finite  element  methods  to  discretize 
moment  models  in  semiconductor  device  simulations.  We  consider  drift-diffusion  (DD)  and  high-field  (HF)  models  of  one 
dimensional  devices,  which  involve  not  only  first  derivative  convection  terms  but  also  second  derivative  diffusion  terms,  as  well 
as  a  coupled  Poisson  potential  equation.  Error  estimates  are  obtained  for  both  models  with  smooth 

solutions.  The  main  technical  difficulties  in  the  analysis  include  the  treatment  of  the  inter-element  jump  terms  which  arise  from 
the  discontinuous  nature  of  the  numerical  method,  the  nonlinearity,  and  the  coupling  of  the  models.  A  simulation  is  also 
performed  to  validate  the  analysis. 


We  have  developed  a  new  cell-centered  control  volume  Lagrangian  scheme  for  solving  Euler  equations  of  compressible  gas 
dynamics  in  cylindrical  coordinates.  The  scheme  is  designed  to  be  able  to  preserve  one-dimensional  spherical  symmetry  in  a 
two-dimensional  cylindrical  geometry  when  computed  on  an  equal-angle-zoned  initial  grid.  Unlike  many  previous  area  weighted 
schemes  that  possess  the  spherical  symmetry  property,  our  scheme  is  discretized  on  the  true  volume  and  it  can  preserve  the 
conservation  property  for  all  the  conserved  variables  including  density,  momentum  and  total  energy.  Several  two  dimensional 
numerical  examples  in  cylindrical  coordinates  are  presented  to  demonstrate  the  performance  of  the  scheme  in  terms  of 
symmetry,  accuracy  and  non-oscillatory  properties. 

We  have  presented  a  high-order  weighted  essentially  non-oscillatory  (WENO)  scheme,  coupled  with  a  high-order  fast  sweeping 
method,  for  solving  a  dynamic  continuum  model  for  bi-directional  pedestrian  flows.  The  dynamic  continuum  model  for 
bi-directional  pedestrian  flows  is  composed  of  a  coupled  system  of 
a  conservation  law  and  an  Eikonal  equation.  We  present  the  first-order 

Lax-Friedrichs  difference  scheme  with  first  order  Euler  forward  time  discretization,  the  third  order  WENO  scheme  with  third 
order  total  variation  diminishing  (TVD)  Runge-Kutta  time  discretization,  and  the  fast  sweeping  method,  and  demonstrate  how  to 
apply  them  to  the  model  under  study.  We  present  a  comparison  of  the  numerical  results  of  the  model  from  the  first  order  and 
high-order  methods,  and  conclude  that  the  high-order  methods  are  more  efficient  than  the  first  order  one,  and  they  both 
converge  to  the  same  solution  of  the  physical  model. 
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